(ns thi.ng.geom.mesh.subdivision
  (:require
   [thi.ng.math.core :as m]
   [thi.ng.geom.core :as g]
   [thi.ng.geom.utils :as gu]
   [thi.ng.geom.vector :as v :refer [vec3]]
   [thi.ng.geom.gmesh :as gm]
   [thi.ng.geom.triangle :as t]
   [thi.ng.dstruct.core :as d]
   [thi.ng.xerror.core :as err]
   [clojure.core.reducers :as r]))

;; Shared helpers

(defn fold-into-map
  [coll] (r/fold (r/monoid into hash-map) conj coll))

(defn face-loop-triples
  "Takes a mesh face (vector of points) and returns lazyseq of
  successive point triples: [prev curr next]"
  [f]
  (->> f (first) (conj f) (cons (peek f)) (d/successive-nth 3)))

;; All subdivision methods as implemented here will currently only
;; function with the default mesh type (GMesh). Support for other mesh
;; types is planned, pending some further refactoring of existing mesh
;; graph functions.
;;
;; Catmull-Clark
;;
;; Based on reference:
;; https://en.wikipedia.org/wiki/Catmull%25E2%2580%2593Clark_subdivision_surface

(defn cc-face-points
  "Takes a seq of faces and returns map with faces as keys and their
    centroids as values."
  [faces]
  (->> faces
       (r/map (fn [f] [f (gu/centroid f)]))
       (fold-into-map)))

(defn cc-subdiv-face
  [face fp e-points]
  (->> (face-loop-triples face)
       (r/map (fn [[p c n]] [(e-points #{p c}) c (e-points #{c n}) fp]))
       (r/foldcat)))

(defn cc-edge-points
  [mesh edges f-points]
  (->> edges
       (r/map
        (fn [e e-faces]
          [e (-> (mapv #(f-points (g/vertices % mesh)) e-faces)
                 (conj (first e))
                 (conj (second e))
                 (gu/centroid))]))
       (fold-into-map)))

(defn cc-subdiv-faces
  [f-points e-points]
  (->> f-points
       (r/mapcat (fn [f fp] (cc-subdiv-face f fp e-points)))
       (r/foldcat)))

(defn cc-replace-vertices
  [mesh f-points e-points sd-faces]
  (let [new-verts (->> mesh
                       :vertices
                       (keys)
                       (r/map
                        (fn [v]
                          (let [f  (->> (gm/vertex-faces* mesh v)
                                        (mapv #(f-points (g/vertices % mesh)))
                                        gu/centroid)
                                vn (gm/vertex-neighbors* mesh v)
                                n  (count vn)
                                r  (gu/centroid (mapv #(m/mix v %) vn))]
                            [v (m/addm (m/madd r 2.0 f) (m/* v (- n 3)) (/ 1.0 n))])))
                       (fold-into-map))]
    (map (fn [f] [(replace new-verts f)]) sd-faces)))

(defn catmull-clark
  [{:keys [edges] :as mesh}]
  (let [faces    (map #(g/vertices % mesh) (g/faces mesh))
        f-points (cc-face-points faces)
        e-points (cc-edge-points mesh edges f-points)]
    (->> (cc-subdiv-faces f-points e-points)
         (cc-replace-vertices mesh f-points e-points)
         (g/into (g/clear* mesh)))))

;; Doo-Sabin

(defn ds-edge-midpoints
  [edges]
  (->> edges
       (r/map (fn [e] [e (m/mix (first e) (second e))]))
       (fold-into-map)))

(defn ds-face-points*
  [f e-points]
  (let [fp (gu/centroid f)]
    (->> (face-loop-triples f)
         (r/map (fn [[p c n]]
                  [c (-> fp
                         (m/+ (e-points #{p c}) (e-points #{c n}))
                         (m/addm c 0.25))]))
         (fold-into-map))))

(defn ds-face-points
  [{:keys [faces edges] :as m}]
  (let [e-points (ds-edge-midpoints (keys edges))]
    (->> faces
         #?(:cljs (seq))
         (r/map (fn [f] [f (ds-face-points* f e-points)]))
         (fold-into-map))))

(defn ds-fp-faces
  "Construct new faces from face points"
  [f-points]
  (map (fn [[f fps]] (mapv fps f)) f-points))

(defn ds-make-edge-loop
  [edges]
  (loop [acc (transient [(ffirst edges)]), p (peek (first edges)), edges (set (rest edges))]
    (if (seq edges)
      (let [e (some #(if (= p (first %)) %) edges)]
        (if e
          (recur (conj! acc (first e)) (second e) (disj edges e))
          (err/throw! "Mesh not manifold")))
      (persistent! acc))))

;; TODO currently only works for GMesh
(defn ds-dual-faces
  [{:keys [vertices] :as mesh} f-points]
  (->> vertices
       (r/map
        (fn [v v-faces]
          (let [vn (gm/vertex-neighbors* mesh v)]
            (if (< 2 (count vn))
              (->> vn
                   #?(:cljs (seq))
                   (r/map
                    (fn [n]
                      (let [f1 (f-points (some #(if (= n (get % :prev)) (get % :f)) v-faces))
                            f2 (f-points (some #(if (= n (get % :next)) (get % :f)) v-faces))]
                        [(get f1 v v) (get f2 v v)])))
                   (r/foldcat)
                   (ds-make-edge-loop))))))
       (r/filter identity)
       (r/foldcat)))

;; TODO currently only works for GMesh
(defn ds-edge-faces
  [{:keys [vertices edges] :as mesh} f-points]
  (->> edges
       (r/map
        (fn [e e-faces]
          (if (= 2 (count e-faces))
            (let [[a b] (seq e)
                  va (vertices a)
                  f1 (f-points (some #(if (= b (get % :prev)) (get % :f)) va))
                  f2 (f-points (some #(if (= b (get % :next)) (get % :f)) va))]
              [(f1 a) (f1 b) (f2 b) (f2 a)]))))
       (r/filter identity)
       (r/foldcat)))

(defn custom-doo-sabin
  [{f* :faces, e* :edges, v* :vertices, :or {f* true, e* true, v* true}}]
  (fn [mesh]
    (let [f-points (ds-face-points mesh)
          f-faces (ds-fp-faces f-points)
          e-faces (if e* (ds-edge-faces mesh f-points))
          v-faces (if v* (ds-dual-faces mesh f-points))]
      (g/into (g/clear* mesh) (concat (if f* f-faces) e-faces v-faces)))))

(def doo-sabin (custom-doo-sabin nil))

;; Butterfly
;;
;; This subdivision scheme will only work with pre-triangulated meshes.
;; The implementation does not check for this for speed reasons. If
;; you're not sure your mesh contains only triangles, pass it through
;; `g/tessellate` first...
;;
;; Reference:
;; http://graphics.stanford.edu/courses/cs468-10-fall/LectureSlides/10_Subdivision.pdf

;; TODO move into tri or gmesh ns?
(defn other-point-in-opp-tri
  [edges a b f]
  (t/other-point-in-tri (first (disj (edges #{a b}) f)) a b))

(defn bf-edge-points
  [edges omega]
  (let [om2 (* -0.5 omega)]
    (->> edges
         (r/map
          (fn [e e-faces]
            (let [[a b] (seq e)
                  [efa efb] (seq e-faces)
                  e1 (m/mix a b)
                  d3 (t/other-point-in-tri efa a b)
                  d4 (t/other-point-in-tri efb a b)]
              (if (and d3 d4)
                (let [d6 (other-point-in-opp-tri edges a d3 efa)
                      d5 (other-point-in-opp-tri edges b d3 efa)
                      d8 (other-point-in-opp-tri edges a d4 efb)
                      d7 (other-point-in-opp-tri edges b d4 efb)]
                  (if (and d5 d6 d7 d8)
                    [e (m/+ e1
                            (m/addm d3 d4 omega)
                            ;;(m/addm d5 [d6 d7 d8] om2)
                            (g/reduce-vector d5 + (fn [x _] (* om2 x)) [d6 d7 d8]))]
                    [e e1]))
                [e e1]))))
         (fold-into-map))))

(defn bf-faces-from-edge-points
  [e-points [a b c]]
  (let [ea (e-points #{a b}), eb (e-points #{b c}), ec (e-points #{c a})]
    [[a ea ec] [b eb ea] [c ec eb] [ea eb ec]]))

(defn butterfly
  ([mesh] (butterfly mesh (/ 8.0)))
  ([mesh omega]
   (let [e-points (bf-edge-points (get :edges mesh) omega)]
     (->> (g/faces mesh)
          #?(:cljs (seq))
          (r/mapcat
           (fn [f] (bf-faces-from-edge-points e-points f)))
          (g/into (g/clear* mesh))))))
